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ABSTRACT 

Numerical simulations of global three-dimensional (3D), self-gravitating discs with a 
gap opened by an embedded planet are presented. The simulati ons are customis ed to 
examine planetary gap stability. Previous results, obtained by iLin fc Papaloizoul from 
two-dimensional (2D) disc models, are reproduced in 3D. These include (i) the de- 
velopment of vortices associated with local vortensity minima at gap edges and their 
merging on dynamical timescales in weakly self-gravitating discs, (ii) the increased 
number of vortices as the strength of self-gravity is increased and their resisted merg- 
ing, and (hi) suppression of the vortex instability and development of global spiral 
arms associated with local vortensity maxima in massive discs. The vertical structure 
of these disturbances are examined. In terms of the relative density perturbation, the 
vortex disturbance has weak vertical dependence when self-gravity is neglected. Vor- 
tices become more vertically stratified with increasing self-gravity. This effect is seen 
even when the unperturbed region around the planet's orbital radius has a Toomre 
stability parameter ~ 10. The spiral modes display significant vertical structure at 
the gap edge, with the midplane density enhancement being several times larger than 
that near the upper disc boundary. However, for both instabilities the vertical Mach 
number is typically a few per cent, and on average vertical motions near the gap edge 
do not dominate horizontal motions. 
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1 INTRODUCTION 

Astroph ysical discs may develop radial structure for several 
reasons (|Armitagdl2uTlh . It has been suggested that proto- 
planetary discs contain 'dead zones' in which the magneto- 
rotational instability is ine fficient, leading to a reduced ac- 
cretion rate in this region (|Gammielll996h . Matter then ac- 
cumulates at the radial boundary of a dead zone and the 
actively accreting region, leading to a local density bump. 

Structure can also be induced by an external potential 
such as an embedded satel lite. A sufficiently mass ive planet 
opens a gap in the disc (|Lin fc Papaloizoul \l 986). and the 
gap edges involve radial structure with characteristic length- 
scales of the local disc scale- height. 

It is well established that localised r adial structure 
in th i n discs can be d ynam ically unstable ([Lovelace et al.l 
Il999l : iLi et al.l l200(j l200lh . The evolution of radially 
structured discs may then be affected by such instabil- 
ities. More specifically, a necessary condition for insta- 
bility is the existence of an extremum in the ratio of 
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vorticity to surface density, or vortensit}Q- Indeed, in- 
stabilities have been demonstrated explicitly for dead 
zone boundaries (IVarniere fc Ta gger 200(1: iLvra et al.ll2008l . 
2009b; Crespe et al. 2 01lh as we l l as planetary gap edges 
(iKoller et al.ll2003|: ILi et al.ll2005l: Ide Val-Borro et al.ll2007l : 
ILvra et al.ll2009al : \Lm fc Papaloizoul 1 20 id ). 

These studies consider non-self-gravitating or weakly 
self- gravitating discs. In fact, instabili ties in structured as- 
trophy sical discs can be traced back to I Lovelace fc Hohlfeldl 
(197 8|), who con si dered self- gravitating particle discs. 
ISellwood fc Kahnl (|l99lh studied a similar system, 
while self-gravita t i ng ga se ous discs were exami n ed by 
IPapaloizou fc Linl (119891). IPapaloizou fc Savoniiel (|l991l ) 
and iMeschiari fc Laughlinl (2008, who adopted disc profiles 
to mimic planetary gaps ). 

ILin fc Papaloizoul (|2011al lbh explored in more detail the 
role of self-gravity on the stability of gaps self-consistently 
opened by a planet. They performed a series of linear and 
nonlinear calculations for a range of disc masses. They found 

1 This is modified by a factor involving the entropy for non- 
barotropic discs. 
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vortex formation in weakly self- gravitating discs and global 
spiral arms in massive discs, but both are associated with 
the gap edge. 

The above studies have employed the razor-thin or two- 
dimensional (2D) disc approximation. It is natural to ex- 
tend these models to three- dimensions (3D). However, note 
that non-axisymmetric instabilities in pressure-supported 
thick discs, i.e. the Papaloiz ou-Pringle instability (PP I ), was 
originally studied in 3 DjPapaloizou fe Pringle! 1 19841 . Il985l . 
Il987l : iGoldreich et al.l Il986h . The vortex-forming instabil- 
ity mentioned above is essentially the PPI operating in a 
thin disc with the density bump being analogous to a torus. 
Som e early studies of slender to r i also included self-gravity 
(e.g. iGoodman &; Naravanlll988l ; I Christ odoulou &; Naravar] 
Il992h . 

Recently, 3D non- self- gravitating, rotationally- 
supported global discs have been simulated with a 
local density bump, eith er set as an initial condition 
(jMeheut et al.ll20ldJ2012bl Iah, or self- cons istently generated 
by a resistivity jump in magnetic discs ([Lyra &; Mac Low! 
2012). The latter models the dead zone scenario. These 
simulations display vortex formation similar to 2D discs. 
On the other hand, instabilities at planetary gap edges in 
3D self- gravitating discs have not yet been simulated. 

I n this work we extend the 2D self- gravitating disc mod- 
els of llLin fe Papaloizoul (|2011al lbh to 3D. Because the insta- 
bilities are associated with radial structure with comparable 
size to the disc thickness, it is not obvious at first that 2D is a 
good approximation. Thus, our priority in this first stu dy is 
to verify results obtained in lLin &; Papaloizo d feonnl lhh by 
simulating equivalent systems in 3D. We also identify some 
three-dimensional effects that sets the direction for future 
investigations. 

This paper is organised as follows. After reviewing the 
main 2D results in the next subsection, we describe our 3D 
disc-planet models in 32] Numerical methods are stated in 
33] We go through our simulations in 31] with additional re- 
sult analyses presented in 3S] We summarise and conclude in 
3H]with a discussion of several limitations of our simulations. 



1.1 Gap stability in 2D 

For discussion purposes here we consider a barotropic disc 
0. For a radially structured disc the quantity governing sta- 
bility is the vortensity profile 
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(1) 



where k 2 = R~ 3 d(R 4 Q 2 ) /dR is the square of the epicy- 
cle frequency, Q is the disc angular velocity and S its sur- 
face density. R is the cylindrical radius. If the disc is self- 
gravitating, the Toomre parameter Q is also important, 



Q = 



c s k, _ c s ( 2Vir\\ 



1/2 



(2) 



2 Our locally isothermal numerical models are not strictly 
barotropic. However, the instabilities of interest are associated 
with localised structure and we adopt sound-speed profiles that 
vary slowly in these regions. Hence we can consider it to be 
isothermal and barotropic. 



where c s is the sound-speed. Q < 1 signifies local axisym- 
metric gravitational instability (jToomrd[l9fri . We remark 
that some studies of instabilities in structured discs employ 
values of k 2 marginally above zero or even negative (e.g. 
ILi et al.ll2000l l2Q0lh . implying the classic Toomre instabil- 
ity may operate. Of course, if the disc is strictly non-self- 
gravitating and k 2 > 0, then £ can be rescaled so that 
Q ^> 1, giving a self-consistent model. 

The connection between Q and 77 results in the vorten- 
sity profile of a planetary gap to resemble its Toomre param- 
eter profile because of vortensit y generation and destru ction 
across planet-induced shocks (|Lin &; Papaloizoul 120101 ). Ex- 
trema in Q and 77 nearly coincide at the same radius. Fig. [T] 
shows a typical Toomre profile. Here we focus on the outer 
disc where we find instabilities strongest in the numerical 
simulations. 

As mentioned in 3H disc profiles with stationary points 
in 77 (and therefore in Q for planetary gaps) can be dy- 
namically unstable. The gap profile fulfils this require- 
ment and has the follow ing stability properties (taken from 
iLin fc Papaloizoul l2011aH bh: 



(i) In weakly or non-self-gravitating discs, instability is 
associated with the vortensity minimum or min(Q), leading 
to local vortex formation. 

(ii) As the strength of self- gravity is increased, the vor- 
tex mode shifts to higher azimuthal wavenumber m. This is 
partly due to the stabilisation effect of self-gravity on low m 
vortex modes. 

(iii) The timescale for vortex-merging increases with the 
strength of self- gravity. In non-self- gravitating discs, vortices 
merge on dynamical timescales and the result is a single, az- 
imuthally extended vortex. Multi- vortex configurations can 
last much longer with increased disc gravity. Merging even- 
tually takes place but the resulting vortex is azimuthally 
lo calised. In the moderately self-gravitating case discussed 
in lLin &; Papaloizoul (|2011a ). a vortex-pair persists until the 
end of the simulation. 

(iv) The vortex mode is suppressed with sufficiently 
strong self-gravity and replaced by a global spiral instabil- 
ity associated with the local vortensity maxima or max(Q). 
The instability can be physically understood as gravitational 
coupling between the gap edge and the wider disc exterior 
to it. 



We shall confirm the above in 3D disc models, except for 
the azimuthally localised, post-merger vortices in (iii) This 
re quires very long simulatio ns with self- gravity (r 



' 200 orbits 

in lLin fc Papaloizoull2011ah , which are currently impractical 
in 3D. 

These instabilities also affect planetary migration, 
leading to vortex-planet and spira l-planet interactions 
(|Lin fe PaDaloizoull20ldl2011al lbL l2012l ). In order to focus on 
gap stability we will not consider migration (and because of 
resolution limits in a 3D simulation), but we can still mea- 
sure disc-planet torques. In particular, we will confirm that 
spiral modes make the disc-on-planet torques more positive 
with increasing instability strength. 
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Figure 1. Typical structure of the outer edge of a planetary 
gap in terms of the Toomre parameter Q. The horizontal axis is 
the displacement away from the planet in units of its Hill radius 
(denned in gZ3) . 



2 DISC-PLANET MODELS 

We consider a three-dimensional gas disc of mass Md with 
an embedded planet of mass M p , both rotating about a cen- 
tral star of mass M*. To describe the system, we use both 
spherical polar co-ordinates r — (r, 0, 0) and cylindrical co- 
ordinates r = (R, 0, z) centred on the star. The frame is 
non-rotating. The disc is governed by the standard fluid 
equations 



dp ^ 
— + V • 

du 



--Vp- 

p 



2 

c s p. 



(pu) = 0, (3) 

V$ eff (4) 

(5) 

where p is the mass density, u is the velocity field, p is 
the pressure and <E> e ff is an effective potential. Physical vis- 
cosity is not included in this study (but artificial viscosity 
is employed in the numerical simulations to treat shocks). 
The effect of viscosity on vortex modes and spiral modes 
have been investigated p reviously (|de Val-Borro et~aL 2007; 
iLin &; PapaloizouT 2011bh . We do not expect its effect to dif- 
fer in 3D. 

2.1 Equation of state 

We adopt a modified isothermal equation of state (EOS) 
where the sound-speed c s depends on R and the planet po- 
sition if present. Without a planet, we set 



No planet, 



(6) 



where H — hR is the disc scale-height with constant aspect- 
ratio h and Qfc = y GM* / R 3 is the Keplerian frequency. 
When a planet is present, we set 



7/2\2/7 



With planet, 



(7) 



+ 4 



where H p — h p d p , — GM p /d p and d p — yjr — r p 
is the softened distance to the planet at position r p and e p 
is the softening l ength defined later. This EOS is taken from 
IPeplinski et al.l (|2008l ) and is used here to increase the tem- 
perature near to planet in order to reduce mass accumula- 
tion in this region. The dimensionless parameter h p controls 
the temperature increase at r p relative to that given by Ci so - 



2.2 Effective potential 

The effective potential is: 
$eff = + ® P + $ + $2, 
where 

GAL 
<£* = 

r 

is the stellar potential and 

GM P 



(8) 
(9) 

(10) 



is the softened planet potential. $ is the gravitational po- 
tential due to the disc material and is given via the Poisson 
equation 



V 2 $ = 4ivGp. 



(11) 



In Eq. [8] &i is the indirect potential due to the disc and the 
planet, 



Mr) 



I 



Gp(r') / 3 / , GM P 

_ r . r d r _|_ -iL T . . r 



(12) 



The indirect potential accounts for the acceleration of the 
co-ordinate origin relative to the inertial frame. This term 
is included for consistency but is unimportant for the insta- 
bilities of interest. 



2.3 Initial disc 

The physical disc occupies r G [n,r ], £ [# m in,7r — 
#min] and (j> G [0, 2tt]. The vertical domain is such that 
tan(7r/2 — O m i n )/h = uh scale-heights. The density field is 
initialised to 



p(t = 0) =fipo, 



(13) 



where po is the density profile corresponding to a non-self- 
gravitating disc, 



po(R,z) 



2ttH 
x exp 



R + hn 



l_ 



(14) 



with fixed power- law index a — 3/2. /3isa function to ac- 
count for vertical self-gravity, such that the surface density 
of the initial disc is the same as that corresponding to po. 
We calculate ft in Appendix [A] with some approximations. 
Because of this, we always first evolve the disc without a 
planet. 

The constant So is chosen via the Keplerian Toomre 
parameter Q Q at the outer boundary: 



Qo 



E 



i: 



R=r a 

p dz. 



(15) 
(16) 



Note that the integration for surface density E is taken over 
the finite vertical domain being considered. 

The disc is initialised with zero meridional velocity 
(u r = uq =0). The azimuthal velocity is set by centrifugal 
balance with stellar gravity, pressure and self-gravity, but 
for the disc models being considered, which are thin and 
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not very massive, the initial azimuthal velocity is essentially 
Kepler ian. 

Our disc models are labelled by Q ocM; 1 . This gives 
an indication of the strength of self-gravity. Specifically it 
measures gravitational stability against local axisymmetric 
perturbations at the outer disc boundary. All of our discs 
satisfy the Toomre criterion for stability. 

2.4 Planet configuration 

In this work the planet is treated as a fixed external poten- 
tial. Its purpose is to create and maintain a structured disc. 
The planet is held on a circular orbit, r p = (r p ,7r/2,0 p ) 
with 4> p (t) = Qk(r p )t in spherical co-ordinates. The soft- 
ening length of the planet potential is fixed to e p = O.lr^ 
where r h = (q/3) 1/3 r p is the Hill radius and q = M p /M*. 
The EOS parameter is set to h p — 0.5. 



3 NUMERICAL METHOD 

We evolve the disc-planet system using the ZEUS-MP finite - 
difference code in spherical coordinates ([Haves et al.ll2Q06h . 
The computational domain is divided into (N r , Nq, N^) 
zones, logarithmically spaced in radius and uniformly spaced 
in the angular coordinates. We assume symmetry about the 
midplane, so the computational domain only covers the up- 
per plane (z > 0). Hydrodynamic boundary conditions are 
outflow at n, r , reflecting at m in and periodic in <j>. 

ZEUS-MP was chosen for its ability to treat self-gravity 
on a spherical grid with parallelisation. It solves the discre- 
tised Poisson e quation using a conjugate gradient method 
(for details, see lHaves et al.ll2QQ6h . To supply boundary con- 
ditions to the solver, we approximate the boundary p oten- 
tial us ing spherical harmonic expansion as described in Boss 
(1980). The expansion in spherical harmonics Yi m is trun- 
cated at /max, ra ma x- We assume negligible contributions to 
the disc potential beyond the physical disc boundaries. 

3.1 Simulation setup 

Computational units are such that G — M* = 1. The 
radial range of the disc is (n,r ) = (1,25). The verti- 
cal extent is nu — 2 scale- heights. The grid resolution is 
(N r , Ne,N(f,) = (256,32,512). We quote time in units of 
Po = 27r/Qfc(r p ). Between ^ t < 10Po the disc is evolved 
without a planet and (/max, r^max) = (48,0). The planet is 
introduced at t = 10Po and its mass smoothly increased 
from zero to its full value between 10Po ^ t ^ 20Po- For 
t > 10P we set (/ max , ^max ) = (16,10). 

The unstable modes of interest are associated with 
vortensity extrema at gap edges, so these features need to 
be resolved. Numerical diffusivity, e.g. due to low resolu- 
tion or grid choice, may inhibit such instabilities. For exam- 
ple, in their 2D studies o f vortex generation at gap edges, 
Ide Val-Borro et all (|2007h did not find instability in Carte- 
sian simulations. 

Test runs with half the resolution in each co-ordinate 
also produce vortex and spiral modes, but with reduced 
growth rate. While the resolution adopted here can confirm 
these instabilities operate in 3D, it may be inadequate to 
probe secondary instabilities on longer timescales (see 36. 3[) . 



Table 1. Simulation parameters. Q p is the Keplerian Toomre 
parameter evaluated at r p . Case was ran without self-gravity. 



Case 


h 


10 3 q 


Qo, Q P 


M d /M* 


mode 





0.07 


2 


oo 


0.021 


vortex 


1 


0.07 


2 


8.0, 14.8 


0.021 


vortex 


2 


0.07 


2 


4.0, 7.40 


0.042 


vortex 


3 


0.07 


2 


3.0, 5.54 


0.056 


vortex 


4 


0.05 


1 


4.0, 7.39 


0.030 


vortex 


5 


0.05 


1 


3.0, 5.54 


0.040 


vortex 


6 


0.05 


1 


1.7, 3.14 


0.070 


spiral 


7 


0.05 


1 


1.5, 2.77 


0.080 


spiral 



We set m m ax = 10 for the boundary potential as a 
compromise between accounting for the non-axisymmetry 
of unstable modes and computation time. We typically find 
vortex modes with azimuthal wave- number m — 3 — 6. We 
tested a run (case 3 below) with m ma x = 12 and obtained 
similar results. However with m ma x = 0, vortex merging pro- 
ceeds soon after their formation, whereas it is resisted with 
m max = 10 (consistent with high-resolution 2D results). The 
global spiral modes have m — 2, so large m m ax is not crucial. 
Indeed, tests with m ma x = 4 also yield the spiral instability. 



4 RESULTS 

Our simulations are summarised in Table [T] If the stel- 
lar mass is taken to be M* = Mq then q — 10 -3 corre- 
sponds to a Jupiter-mass planet. The planetary masses con- 
sidered here are larger th an our previous 2D investigations 
(|Lin &; Papaloizoull2011 a,b) in order to achieve higher insta- 
bility growth rates and shorten the computation time. We 
will examine gap stability as a function of Qo- 

4.1 Vortex modes in weakly self-gravitating discs 

We first compare Case and Case 1. The setup for these 
runs are identical except that the disc potential is neglected 
in Case 0, which corresponds to the standard approach t o 
model disc-planet systems (e.g. iD'Angelo &; Lubowll201oh . 
Case 1 is the self- gravitating version of Case 0. We show 
below that, even with a minimum Toomre parameter of 
Qo = 8, disc self- gravity affects the evolution of the vor- 
tex instability. 

Fig. [2] shows the midplane gap profile at t = 25Po in 
terms of the relative density perturbation. The snapshot 
is taken before instabilities develop. Case 1 has a slightly 
deeper gap and steeper gap edges than Case 0. This is be- 
cause in a self- gravitating calculation such as Case 1, fluid 
bound to the planet adds to its mass. For both runs, the 
mass in the planet's Hill sphere is ~ 0.02M P . 

4-1.1 Development of instability 

Fig. [3] shows the evolution of the gap. It is clear that the 
vortex instability can develop in 3D with or without self- 
gravity. At t = 30Po , two vortices are visible in Case while 
3 vortices are seen in Case 1 (in both cases there may be 
another vortex coinciding with the outer planetary wake). 
For the weakly self- gravitating discs considered here, 
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Figure 3. Development of the vortex instability in the non-self- gravitating Case (top) and the weakly self- gravitating Case 1 (bottom). 
The midplane relative density perturbation is shown. White lines at intermediate azimuths corresponds to the vortex centroid in Fig. 0] 
while the azimuthal range marked by the lines was set for averaging in Fig. [5] 
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Figure 2. Gap profile opened by a giant planet in a non-self- 
gravitating disc (solid, Case 0) and a self-gravitating disc (dotted, 
Case 1). The azimuthally-averaged relative density perturbation 
in the midplane is shown. Profiles for other cases are similar. 

vortex modes with the same m have been excited, but with- 
out self- gravity vortices merge soon after formation. (In Case 
0, the over-density at the outer gap edge just ahead of the 
planet appears to be a merging vortex-pair, rather than a 
single vortex from the instability.) At t = 40Po, only a 
vortex-pair remains in Case while merging is delayed in 
Case 1 with a 3- vortex configuration. 

The differences between self- gravitating and non-self- 
gravitating simulations shown in Fig. [3] are similar to those 
seen in 2D simulations (jLvra et al. 2009b; iLin &; Papaloizoul 



l2011aT ). ILin &; Papaloizoul demonstrated that self-gravity 
sets a minimal inter-vortex distance so that merging is re- 
sisted. Three-dimensionality does not affect the evolution of 
the vortex instability in the rcj) plane. However, at this level 
of self- gravity merging is only delayed. The final state for 
both cases is a single vortex circulating the gap edge. 

4-1.2 Effect of self- gravity on vortex vertical structure 

Here we compare the final vortex in Case and Case 1. Fig. [4] 
shows their vertical structure in the Rz plane. The snapshots 
correspond to the vortex centroid, marked by white lines at 
intermediate azimuths in the right panel of Fig. [3] Without 
self-gravity the vortex instability produces predominantly 
columnar disturbances in the relative density perturbation 
W. 

Case is consistent with recent linear calculations of 
the vortex instability in non- self- gravitating 3D discs, which 
show that W has essentially no vertical depe ndence at the 
radius where vo rtex-formation is expected (Me heut et al.l 
l2012bl : lLm1l2012h . In Fig. g] the Case vortex does show 
very weak vertical dependence near the upper boundary. 
This is likely due to the finite vertical domain adopted in 
our modejf]. 

3 Linear calculations of vertically isothermal discs usually assume 
an atmosphere of infinite extent. 
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Figure 4. Vertical structure of a vortex in a non-self-gravitating 
disc (top, Case 0) and a weakly self-gravitating disc (bottom, 
Case 1). Contours show the relative density perturbation and the 
arrows are mass flux vectors pu projected onto this plane. The 
snapshots correspond to the vortex centroids marked by white 
lines at intermediate azimuths in Fig. [3] (right panel). 

By contrast, the self- gravitating vortex in Case 1 clearly 
display stratification in the relative density perturbation. 
The vortex is more concentrated toward the midplane. With 
self-gravity, density enhancement in the vortex can be ~ 
50% higher than without. 

In Fig. \E\ we plot vertical velocities averaged over the 
vortices. The non- self- gravitating case typically involve s pos- 
itive vertical velocity (also seen in Me heut et al ]l2012bh . Per- 
haps not surprisingly, the self- gravitating case has as strong 
over-density near the midplane to provide vertical acceler- 
ation, so on average the vertical velocity is negative. The 
precise values in Fig. \E\ depends on the averaging procedure 
but the contrast in sgn(i^) between the two cases is robust 
(even when we consider the point in the Rcj) plane where den- 
sity perturbation is largest and do not perform an average). 
The quantity (u 2 z / 'cf so ) 1//2 behave similarly. 

Although there is vertical motion, it is worth noting 
that the vertical Mach number is only a few per cent. Fig. 
[3] also indicate that the flow is horizontal on average. This 
suggests approximate vertical hydrostatic balance. If self- 
gravity is included, then just like in the set up of the initial 
disc, the additional vertical force will enhance the midplane 
density. Hence, we observe a more stratified vortex in Case 
1. 

4.2 Vortex modes with moderate self-gravity 

Cases 2 — 5 are all self- gravitating and develop the vortex 
instability. Case 2 and Case 3 are continuations of Case 1 
to more massive discs. Case 4 and Case 5 are identical runs 
except a Jupiter-mass planet (q = 10 -3 ) and a thinner disc 
(h = 0.05) are adopted. 

We compare the relative density perturbation between 
Case 2 and Case 3 in Fig. \6\ The snapshots are taken at 
z — H but are very similar to razor-thin disc simulations 
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Figure 5. Average vertical velocity u z , normalised by the sound 
speed, in a non-self-gravitating vortex (Case 0, solid) and a self- 
gravitating vortex (Case 1, dotted). The radial range for the aver- 
age is taken over r — r p £ [4, 6.5]r^ for Case and r — r p G [4, Gjr^ 
for Case 1 since the latter vortex is slightly smaller (see Fig. 3]). 
The azimuthal range is that marked by white lines in Fig. [3] (right 
panel). 

(|Li et al.ll2009l : lYu et al.lboiol : iLin fe Papaloizoull2011ah . In 
both Cases the m — 5 vortex mode is excited (cf. m — 3 — 
4 in Case 1). At t = 30Po, vortices are just beginning to 
emerge in Case 3, whereas distinct blobs can already be 
identified in Case 2, with larger over-densities than vortices 
in Case 3. This indicates a stronger instability with decreas- 
ing strength of self-gravity. Vortex merging ensues in Case 
2 but does not occur in Case 3 within the simulation time- 
scale (unlike Case 1). 

For razor-thin discs, ILin &; Papaloizoul (|2011ah have 
shown that self-gravity has a stabilising effect against vortex 
modes with low m. This contributes to favouring higher m 
modes and hence more vortices with increasing self-gravity. 
Fig. [6] together with results for Case 1, show that this effect 
persist in 3D. Resisted-merging, seen in 2D models, also oc- 
cur in 3D. Because this is due to vortices executing mutual 
horseshoe turns, only horizontal self- gravitational forces are 
important. We do not expect the vertical dimension to sig- 
nificantly affect vortex- vortex gravitational interaction. 



4.2.1 A 3D vortex in Case 3 

Here we examine the vortex in Case 3 marked by horizontal 
white lines in Fig. [6] Note that no merging has occurred. 
The vortex is radially located about r v ~ r p + 5r^. It has 
azimuthal and radial sizes A(f) v ~ 9h and Ar v ~ 2.5if , 
respectively. Its mass is M v ~ 8.4 x 10- 4 MH The vortex 
Hill radius ru v ~ 0.9H is smaller than its horizontal size but 
comparable to its vertical size at the vortex centroid. 

Fig. [3 shows the vertical structure of the vortex de- 
scribed above. As expected it is more stratified than in 
weakly self- gravitating discs (Fig. [4}, especially at the vor- 
tex centroid where the over-density is maximum. The initial 
Kepler ian Toomre parameter at r = r v is Q = 4.3. A density 



4 Since no merging has occurred, M v can also be estimated by 
the mass contained in an annulus about r v in the unperturbed 
disc divided by m = 5. This gives 8 x 10 -4 M* if the annulus 
width is 2r/j,. 
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Figure 6. Development of the vortex instability at gap edges in 
moderately self- gravitating discs. The relative density perturba- 
tion is shown for Case 2 (top) and Case 3 (bottom) at z = H . 
The planet position is marked by the cross. Horizontal lines in 
the plot for Qo = 3, t = 50Pq indicate azimuths taken in Fig. 



enhancement by a factor > 2 easily gives Q < 2, so that 
self- gravity in the perturbed state is dynamically important. 

Ahead and behind the vortex centroid, the flow field is 
mostly horizontal (top and bottom panels in Fig. [7]). It corre- 
sponds to the moti on of an ant i- c yclonic disc vortex common 
in 2D simulations (jLi et al.ll200ll ). In this plane, vertical mo- 
tions only become significant compared to ur close to the 
vortex centroid (middle panel). We plot the average u z for 
the vortex in Fig. [8] Azimuthal slices at the vortex centroid, 
behind it and ahead of it are also shown. The average vertical 
velocity is positive (solid line). This is qualitatively differ- 
ent to the merged vortex in the weakly self- gravitating Case 
1 (Fig. 2]). We have examined other vortices but could not 
identify a 'typical' vertical flow structure (c.f. anti-cyclonic 
horizontal flow is generic). Although the vertical velocities 
appear to display a range of behaviour, the vertical Mach 
number is still very small. 

We also find the vortices have height-dependent az- 
imuthal structure. Fig.[9]shows several slices in the <j>z plane. 
The choice of radii for these plots were based on the vortex 
described above (corresponding to the right vortex in the fig- 
ure) . The vortices are more columnar closer to the gap edge 
(r — r p = 4rn, top slice). Moving away from the gap edge, 
the vortices develop front-back asymmetry and are thinner 
with increasing height (r — r p = 5.5r^). 

The increased three-dimensional structure away from 
the gap edge could be related to dens ity waves emitted 
by the vortex ([Paardekooper et al.l 12010). Back-reaction of 
these waves on the vortex vertical structure, if any, is ex- 
pected to be weaker on the side of the vortex adjacent to 
the low-density gap edge (top slice in Fig. [9}. 

The perturbed azimuthal velocities away from the vor- 
tex centroids again follow the expected anti-cyclonic motion 
in the horizontal plane (being positive interior and nega- 
tive exterior to the centroid). However, unlike the previous 
Rz plots, here the vertical velocities can often be compa- 
rable or larger than the perturbed azimuthal velocities (e.g. 
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Figure 7. Vertical structure of a vortex in the moderately self- 
gravitating Case 3. The relative density perturbation are overlaid 
by the mass flux vectors pu projected onto this plane. The slices 
are taken at t = 50Po and azimuths ((f) — (j) p ) / ir = 0.86 (top), 0.81 
(middle) and 0.76 (bottom). These correspond to the horizontal 
lines marked in Fig. [6] 
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Figure 8. Vertical velocities u z , normalised by the sound speed, 
for the vortex in the self-gravitating Case 3 shown in Fig. [6] The 
velocities are all radially averaged over r — r p G [4.0, 5.5]r^.The 
solid line is also averaged over the vortex azimuthal extent. 
Azimuthal slices ahead of the vortex (dotted), at the centroid 
(dashed) and behind the vortex (dashed-dot) are also shown. 



r — r p — 5.5r^, bottom panel). Of course, the total azimuthal 
velocity is supersonic and therefore much larger than verti- 
cal velocities. As remarked above, the vortices do not share 
the same structure. For example, the left vortex centroid in- 
volves negative vertical velocity while u z > for the right 
vortex (r — r p = 4.75r^, middle panel). 
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Figure 9. Azimuthal vortex structures in the self-gravitating 
Case 3 at t = 50Po. The background flow is from left to right. 
The relative density perturbation are overlaid by the perturbed 
mass flux vectors p(u — RQk<fi) projected onto this plane. The 
slices are taken at radii (r — r p ) /r^ =4.0 (top), 4.75 (middle) 
and 5.5 (bottom). These correspond to the flow just interior, at, 
and just exterior to the vortex centroid on the right. 



4-2.2 Case 4 and Case 5 

Cases 4 — 5 are additional examples of the vortex instabil- 
ity in 3D discs with smaller h (= 0.05) than the above runs 
(with h — 0.07). However, initial Toomre profile is nearly in- 
dependent of h, so the strength of self- gravity is unchanged. 
We also used a smaller planetary mass, q = 10" 3 so the 
ratio q/h 3 = 8 is not significantly larger than Cases 2 — 3 
(q/h 3 — 5.8). We checked that prior to instability, the outer 
gap edge profiles of these cases are similar. 

We found the instability grows slower with h = 0.05 
as vortices become identifiable at t = 35Po, compared to 
t = 30 Po for Cases 2 and 3. A decrease in growth rate with 
sound-sp eed (which is proportional to h) was already noted 
m 2D (|Li et al. 2000). Thus, the effect of sound-speed on the 
vortex instability remain unchanged by the 3D geometry. 

Fig. [10] show several snapshots of Case 4 and Case 5 
at t — 50Po near the upper disc boundary. As before, high 
in the atmosphere the density enhancement is weaker with 
increasing self- gravity. Consistent with 2D simulations, the 
more self- gravitating Case 5 develops the m — 6 vortex mode 
whereas m = 5 vortices develop in Case 4 with weaker self- 
gravity. Merging is strongly resisted in these runs. In fact, 
Case 5 was extended to t = 135Po and only one vortex-pair 
merged. 
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Figure 10. Vortex instability in Case 4 (left) and Case 5 (right). 
The relative density perturbation near the upper vertical bound- 
ary . These runs are similar to Cases 2 — 3, except a colder disc 
is employed (h = 0.05) and a smaller planetary mass is used to 
open the gap q = 10 -3 . 

4.3 Edge-spiral modes in massive discs 

iLin fc Papaloizoul (|2011al lbt) found that as the strength of 
self- gravity is increased, instability eventually shifts from 
localised vortices to low m global spirals extending from the 
outer gap edge to the outer disc. The vortex instability is 
suppressed because they are asso ciated with vortens i ty min - 
ima (which we will check in flOl . lLin fe Papaloizoul (|2011aT ) 
gives a simple energy argument which show that such as- 
sociation is not possible for sufficiently strong self-gravity. 
Association with vortensity maxim um is favoured in stead. 
These are the spiral instabilities. ILin &; Papaloizoul called 
them edge modes since they can still be considered associ- 
ated with the gap edge. 

We begin to observe edge modes in our 3D models when 
Qo = 1.7 (Case 6)0. Its evolution is depicted in Fig. [TT] 
Am — 2 — 3 disturbance develops at the outer gap edge 
at t = 35Po and induces spiral density waves in the exte- 
rior disc through self-gravity. Interaction between the edge 
disturbance and the wider disc leads to the global spiral 
pattern seen at t = 40Po and t = 45Po. This coupling is 
necessary for instability. The edge mode is associated with 
a local max(Q) just inside the unperturbed outer gap edge. 
In Case 6, max(Q) ~8soa local gravitational instability is 
not possible. However, we can still consider the global edge 
mode as being composed of two parts: an edge disturbance 
where density perturbation is largest, and the spiral arm it 
induces. 

Comparison between the two heights in Fig. [TT] show 
that edge modes are significantly vertically stratified. Most 
of the density perturbation is confined near the midplane. 
Unlike vortex modes the spirals appear transient. Their am- 
plitudes are much reduced by t = 50Po, but are still visi- 
ble. This is likely a radial boundary effect. The transition 
to a low density annulus towards the outer disc edge is a 
result of the standard outflow boundary condition applied 
there. While the linear edge mode instability is insensitive 
to boundary conditions, its long term evolution is affected. 

5 Note that the transition from vortex to spiral modes with de- 
creasing Qo is not abrupt. It is possible to have a mixture. 
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Figure 11. Case 6: development of the edge mode instability at the gap edge of a Jovian mass planet. The relative density perturbation 
in the midplane (bottom) and in the atmosphere (top) is shown at times (left to right) t = 35Po, 40Po, 45Po, 50Po. White lines indicate 
azimuthal cuts taken in Fig. 1121 



The outer disc edge can reflect waves back to to gap edge to 
stabilise it, causing saturation (|Lin &; Papaloizou1l2011bl ). 



4-3.1 Vertical structure of an edge mode 

Several vertical cuts of the edge mode in Case 6 are shown in 
Fig. 1121 The slices are taken at azimuths marked by white 
lines in Fig. [11] The top three plots are associated with 
the edge disturbance, while the bottom plot is taken at the 
transition between the edge disturbance and its spiral arm 
extending to the outer disc. 

Fig. [12] shows that the horizontal flow in the edge dis- 
turbance differ significantly from the vortex mode. If the 
second plot is considered the 'centroid', then the inward 
(outward) flow ahead (behind) it is in the opposite sense 
to anti-cyclonic motion associated with a vortex mode. The 
centroid is the most stratified region. Its has midplane over- 
density rsj 3.37 corresponds to a Toomre Q ~ 1 in the cen- 
troid, but we do not observe fragmentation. 

As we move away from the centroid in the decreasing 
<f> direction, the edge mode decreases in amplitude but oc- 
cupies more of the vertical domain. The bottom plot in 
Fig [12] shows a radial split: the columnar disturbance in 
R — r p G [5.8, 6.5] ru is the beginning of the spiral wave ex- 
cited by the edge disturbance. The spiral arm is fully three- 
dimensional. 

In Fig. [13] we plot the average vertical velocity inside 
an edge disturbance of the spiral mode. A simulation with 
Qo = 1.5 (Case 7) is also plotted for comparison (its mid- 
plane density perturbation is shown in Fig. [T4j) . The flow is 
typically downwards toward the midplane. This is expected 
because of strong vertical self-gravity, as these are massive 
discs and the midplane density enhancement is large. Note 
that u z approaches zero again beyond z/H ~ 1.5 because of 
the imposed reflective upper boundary. 



5 ADDITIONAL RESULTS ANALYSIS 

In this section we examine some secondary quantities de- 
rived from the hydrodynamic simulations above. To keep 
this discussion concise, we will use selected simulations from 
above for illustration. 



5.1 Three-dimensionality 

A simple measure of three-dimensionality of the flow is to 
compare vertical to horizontal motion. Since we are inter- 
ested in non-axisymmetric perturbations to the gap edge, 
we first Fourier transform the meridional momentum densi- 
ties 



p2tt 

(vRm,v zm )= / p x (uR,u z )exp(-im(j))d(/). (17) 
J 

We define the three dimensionality as ®m(z/H), where 

2 (\v zm \ 2 ) 



e. 



(\v zm \ 2 ) + (\v R m\ 2 y 



(18) 



and (•) denotes a radial average. Admittedly, this is a crude 
measure, and exact values of G m varies somewhat with de- 
tails of the average. However, we have experimented with 
different averaging domains and found the features described 
below are robust. 

The top panel in Fig. [15] shows G m for Cases 1-3 at 
t = 40Po- The radial average is taken over r — r p £ [3, 7]rn- 
These are all vortex modes (see Fig. [3]and Fig.©. The flow 
becomes increasingly three-dimensional away from the mid- 
plane but G m = O(10 _1 ) is small. In an averaged sense the 
flow is mostly horizontal. At the end of the simulation for 
Case 1, an azimuthally extended vortex dominates the flow, 
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Figure 12. Vertical structure of the spiral instability at the gap 
edge in Case 6. The relative density perturbation are overlaid by 
the mass flux vectors pu projected onto this plane. The slices are 
taken at t = 40Po and azimuths marked by white lines in Fig. 1111 
(with decreasing <j> from top to bottom). 




Figure 13. Average vertical velocity of the edge disturbance in 
Case 6 (solid) and Case 7 (dotted). The average is taken over 
r — r p £ [2.5, 5.5]r^. The azimuthal range is indicated by white 
lines in Fig. [l2]for Case 6 (first and third azimuth in the clockwise 
direction) and in Fig. [14] for Case 7. 
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Figure 14. The m = 2 edge mode in Case 7. The colourbar range 
is the same as in Fig. 1111 White lines indicate the azimuthal 
range taken for averaging in Fig. [13] (dotted line). The relative 
density perturbation of the edge disturbance in the Rz plane is 
very similar to Case 6 (Tig. I12|) . 



for which we measured Bi ~ 0.2 — 0.3. Thus, although verti- 
cal motion can become an appreciable fraction of horizontal 
motion, the former never dominates. 

G m for Cases 4 — 7 are shown in in bottom panel of Fig. 
1151 The radial average is performed over r — r p £ [2, 6]r^ 
because the global spirals in Cases 6 — 7 significantly pro- 
trude the gap edge. The snapshot is taken at t = 50Po for 
Cases 4 — 5, at t = 40Po for Case 6 and at t = 30Po, so 
that the vortices and spirals have comparable over-densities 
at the gap edge. It also reflects the fact that spiral modes 
are more unstable than v ortex modes and develop earlier 
(|Lin &; Papaloizoul l2011bh . G m ~ 0.2 is again not particu- 
larly large, but the spiral modes are distinctly more three- 
dimensional than vortex modes. This is likely due to ad- 
ditional vertical acceleration provided by the strong self- 
gravity in those cases. 

5.2 Vortensity field 

A fundamental distinction between the linear vortex and 
edge mode instability is their association with local vorten- 
sity minimum and maximum, respectively. In this section we 
compare vortensity fields of discs with vortex modes (Case 
2) and edge modes (Case 7). More specifically, we exam- 
ine the relative perturbation to the vertical component of 
vortensity, 

rj z (t = 0) 



Ar) z 



Vz(t = 0) 



where 



z • V x u 



(19) 



(20) 



Fig. [16] compares Ar] z in the midplane when vortices 
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Figure 15. Three-dimensionality of the non-axisymmetric flow 
near the outer gap edge. Top: Cases 1 — 3 (vortex modes). Bot- 
tom: Cases 4 — 5 (vortex modes) and Cases 6 — 7 (spiral modes). 
The azimuthal wavenumber m is chosen to match the number of 
vortices or large-scale spirals observed. 



and spirals develop. For planetary gaps, vortensity maxima 
and minima are both located near the gap edges with charac- 
teristic separation of the local scale- height. The vortensity 
ring at the inner gap edge (r — r p c± —2rh) remain well- 
defined. The vortex instability is associated with the local 
vortensity minimum near the outer gap edge — seen as lo- 
calised closed contour lines centred about r — r p ~ Arn- The 
vortensity ring at r — r p ~ +2r^ becomes distorted as a con- 
sequence of large-scale vortex formation just exterior to it. 
By contrast, the edge-spiral mode is associated with the lo- 
cal vortensity maximum. Their development inherently dis- 
rupts the vortensity rings. This is seen in the right panel as 
the outer ring is broken up. 

The vertical structures also differ. Fig. 1 1 71 compares Arj z 
at azimuths coinciding with a vortex or the edge disturbance 
of the spiral mode. Both instabilities involve Arj z < 0. It is 
clear that the spiral mode has stronger vertical dependence. 
Its region of Ar] z < becomes thinner away from the mid- 
plane. In the vortex case this region remains about the same 
width and Ar] z is approximately uniform within it. While 
min(A77 2 ) is of comparable magnitude, the vortensity ring 
at r — r p = 2r^is much weaker and thinner in the spiral case 
(Ari z being a factor ~ 4 smaller than the vortex case). 



5.3 Disc-planet torques 

The presence of non-axisymmetric disturbances at gap edges 
is expected to significantly affect disc-planet torques. It has 
been confirmed in 2D simulations that both vortex and spi- 
ral modes lead to oscillatory to rques of either sign (|Li et al.l 
2005; [Lin &; Papaloizoull2011bh . It this section we measure 




Figure 16. Relative perturbation to the vertical component of 
vortensity at the midplane in a disc with the vortex instability 
(left, Case 2 at t = 40Po) and the spiral instability (right, Case 7 
at t = 30Po)- Negative perturbations in the region r — r p G [2, 6]rh 
are outlined by white lines. Dotted horizontal lines indicate az- 
imuthal cuts taken in Fig. ITT1 
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Figure 17. Relative vertical vortensity perturbation associated 
with the vortex instability (top) and spiral instability (bottom). 
The slices are taken at azimuths shown by white dotted lines in 
Fig- EH Negative perturbations in the region r — r p G [2,5.5]r^ 
are outlined by white lines. 



the disc-on-planet torques in several of the above simulations 
to confirm the main features found in 2D. 

We calculate the specific torque acting on the planet 
due to a mass element as 



dT(r) 
f(r,r p ) 



r p x rGp(r)d 3 r 



dl 



exp 



e c r h 



(21) 



(22) 



The tapering function / reduce contributions from close to 
the planet, thereby reducing numerical artifacts arising from 
this region because of the diverging potential and limited 
resolution. We set the parameter e c = 1 so that tapering 
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Figure 18. Instantaneous disc-on-planet torques in simulations 
where the vortex mode develops. Top: Case 3. Bottom: Case 5. 
Note that Case 5 has been extended to t = 135Pq. 



does not significantly reduce contributions from the insta- 
bilities, since they develop at > 2rn away from the planet's 
orbital radius. 

Fig. [18] shows the disc-on-planet torques in Case 3 and 
Case 5, which develop the m — 5 and m — 6 vortex modes, 
respectively. These plots ar e qualitatively similar to 2D sim- 
ulations fe.g. lLi et al.ll2005h . The torques oscillate on orbital 
time-scales and its instantaneous values can be of either sign. 
However, upon averaging over the simulation we find the to- 
tal torques are negative in both cases. This means inwards 
migration is still favoured. 

We extended Case 5 to t = 135Po and find the vortices 
have similar over-densities as at t = 50Po- However, Fig. [18] 
show the torque oscillation amplitudes decrease towards the 
end of Case 5 compared to t G [40, 80]P . At t = 50P the 
vortices are located in r — r p £ [3.5, 5.5] rn but by t = 135Po 
they are located in r — r p G [4, 6]r^. Given that t G [40, 80]Po 
is only 20Po to 50Po after the planet potential has been fully 
introduced, gap- formation is probably ongoing during this 
time. We expect torque amplitudes to be larger during gap- 
formation since the vortices lie closer to the planet. 

Next we examine disc-planet torques in the presence of 
the spiral modes. The top panel of Fig. [19] shows the in- 
stantaneous disc-on-planet torques. Contributions from the 
inner disc (r < r p ) and outer disc (r > r v ) are plotted sep- 
arately for comparison with Fig. 18b in [Lin &; Papaloizoul 
(|2011bh , which is similar to the present plot. Large oscilla- 
tions in the outer torque due to edge mode spirals cause the 
total torque to be positive or negative at a given instant. 

Unlike the vortex modes, Fig. [19] shows that spiral 
modes can lead to a positive running-time averaged torques 
(bottom panel). The average torques become more posi- 
tive with time after spiral modes develop, and with increas- 
ing self-gravity (which increases the instability strength). 
This was also observed in high-resolution 2D simulations in 
iLin fc PapaloTz^ul (|2011bl ). There it was suggested that the 
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Figure 19. Disc-on-planet torques in the presence of spiral modes 
associated with the outer gap edge. Top: total torque (solid), 
torque from the inner disc (dotted) and from the outer disc 
(dashed) in Case 7. Bottom: time-averaged torques in Case 6 
(solid) and Case 7 (dotted). 



creation of large 'voids' in between spiral arms decreases the 
time-averaged density in the planet-induced wakes, thereby 
reducing associated torque magnitudes. Since the outer 
planetary wake normally provide a negative torque, the spi- 
ral modes make the total torque more positive. The sim- 
ilarity between 2D and 3D results indicate that outwards 
migr ation induced by spiral mode s, which was seen in 2D 
bv lLin fe Papaloizoul (|2011bl . 120121 ) , will also operate in 3D. 



6 SUMMARY AND DISCUSSION 

We have performed customised numerical simulations of 
three-dimensional self- gravitating discs, in which an em- 
bedded satellite or planet has opened a gap. We explicitly 
verified in 3D the main results on gap stability previously 
obtained from 2D calculations (iKoller et all 20031; Li et all 
l2005l : lde Val-Borro et al.ll2007l: iMeschiari fe Laughlinll2008l ) . 
in particular those bv lLin Papaloizoul (|2011al lb[). 

Planetary gaps are potentially unstable because of the 
existence of vortensity extrema generated by planet-induced 
shocks. Disc-satellite interaction also occ ur in other systems 
such as stars in black hole accretion discs ([Kocsis et al.ll201ll ; 
iMcKernan et al.ll201ll : iBaruteau et al.ll201ll ~ Furthermore, 
gaps opened by satellites are just one type of structured disc. 
Other examples include dead zone bound aries mentioned in 
JT] and transition discs (Reg alv et al.ll20l3 ). 

Thus, although we were motivated by previous works 
on planetary gap stability, and hence considered disc-planet 
systems, we expect our results to be applicable to discs with 
radial structure of other origin, provided the vortensity pro- 
files involve stationary points and therefore prone to the 
same instabilities. 
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6.1 Confirmed 2D results 

We demonstrated the development of the vortex instabil- 
ity at the outer gap edge opened by a giant planet in 3D 
discs. We began with a non-self-gravitating disc, in which a 
few vortices develop then quickly merge. The quasi-steady 
state is a single azimuthally extended vortex at the gap 
edge. This evolution is similar to the 2D simulation in 
Ide Val-Borro et al.l (|2007T ). 

We also showed the effect of self- gravity on the 
vortex instability obse rved in 2D (iLvra et al. 2009b; 
iLin &; Papaloizoul l2011ah persists in 3D. We observe more 
vortices as the strength of self-gravity is increased by in- 
creasing the density scale. These vortices resist merging on 
dynamical timescales. In our disc models with Qo = 3, the 
multi- vortex configuration lasts until the end of the simula- 
tions. 

As expected from 2D linear theory (|Lin &; Papaloizoul 

l2011bh . vortex modes are suppressed in our massive disc 
models. Instead, a global spiral instability develops which 
is associated with the local vortensity maximum just inside 
the outer gap edge. These are distinct from vortex modes 
since self- gravity is essential. 

Our limited numerical resolution does not permit ac- 
curate disc-planet torque measurements, but the qualitative 
effect of the vortex and spiral instabilities, previously stud- 
ied in 2D, have been reproduced in 3D — oscillatory torques 
of either sign and the tendency for spiral modes to provide 
on average a positive torque. The similarity to 2D results is 
not surprising since for giant planets ru ^ H, so the razor- 
thin disc approximation is expected to be valid as far as 
disc-planet interaction is concerned. Furthermore, vertical 
self-gravity increases the midplane density while reducing 
that in the atmosphere (Appendix |A|), so that given a fixed 
temperature profile the 2D approximation is even better for 
self- gravitating discs. 

It is worth mentioning here that i n shearing sheet sim- 
ulations, iMamatsashvili &; Ricel (|2009l ) found self-gravity to 
favour vortices of smaller scale. They initialized a local patch 
of an unstructured disc with random velocity perturbations. 
Their gravito-turbulent discs are dominated by small vor- 
tices limited by the local Jeans length (which is smaller than 
the scale- height). Without self-gravity, they merge to form 
larger vortices. These observations are similar to the above 
results for the vortex mode. However, the setups are quite 
different as we consider radially structured, laminar global 
discs. Our large-scale vortices develop from a linear insta- 
bility and have horizontal sizes of a few scale-heights. Thus, 
confirmation of the above results in 3D is only valid for the 
edge instabilities considered in this paper. 

6.2 Effects of 3D 

In our simulations the dominant three-dimensional effect is 
vertical self- gravity on the density field. In the non-self- 
gravitating limit, the relative density perturbation associ- 
ated with a vortex is columnar with weak vertical depen- 
dence. This is consistent with w i th 3D lin ear and nonlinear 
simulations (|Meheut et alJl20lA l2012bl lal: lLmll2012h . 

As the strength of self- gravity is increased, vortices be- 
come more vertically stratified — they are condensed to- 
wards the midplane. For moderately self- gravitating discs, 



the vortex midplane density enhancement can be twice that 
near the upper boundary. The spiral modes display signifi- 
cant vertical structure near the gap edge, while the density 
waves they launch in the outer disc is columnar. The latter 
is probably due to the chosen equation of state (see below). 

The effect of vertical self-gravity on the vortex mode is 
seen even in our least massive disc model with Qo = 8. One 
can consider an initially smooth disc which is justified to be 
non- self- gravitating. However, this approximation may be- 
come less good with the creation of a vortensity minimum, 
because it can also be a local minimum in the Toomre Q 
(Eq. [2}. That is, the Toomre parameter is decreased with 
the development of local radial structure (such as a den- 
sity bump). The non- self- gravitating approximation worsens 
further when the vortex instability associated with min(?7) 
develops, because the vortices are regions of enhanced den- 
sity, especially if they merge into a single large vortex. Thus, 
the non- self- gravitating approximation is not guaranteed to 
hold in the perturbed state even if it does in the initial disc. 

It is interesting to note that linear calculations of 
the vortex instability in vertically isothermal , non-self- 
gravitating discs show that t he vertical velocity van- 
ishes near the vortex centroid ([Meheut et al.l l2012bl ; iLinl 
l2012h. but this is not observed in nonlinear calculations 
([Meheut et al.|[2012bL and in the present simulations) . This 
contrasts to their anti-cyclonic horizontal flow, which can be 
compute d in linear theory an d seen in hydrodynamic sim- 
ulations (|Li et alJlioOCl l200lh . This suggests that vertical 
motions in this case may be asso ciated with secondary pro- 
cesses (e.g Goodman et al. 1987|). 

Although we observe somewhat complicated vertical 
flow for both vortex and spiral instabilities, the vertical 
Mach number is at most a few per cent in magnitude. Also, 
the magnitude of vertical motion is at most ~ 20% of the ra- 
dial motion on average. This suggests that the disturbances 
at the gap edge is roughly two-dimensional in the present 
disc models. This would be consistent with early stud- 
ies which find instabilities asso ciated with co-rotation sin- 
gularities are t wo-dimensional (Papaloizou & Pringle 1985; 
iGoldreich et al.l 1 19861 ) . Recent 3D linear calculations also 
find that, even with a vertical temperature gradient, the 
vortex mode (witho ut se lf-gravity) is largely 2D near the 
vortensity minimum (|Linll2012h . 

6.3 Caveats and outlooks 

In order to provide 3D examples of previous 2D results, a 
range of disc models had to be simulated, each for many 
dynamical timescales with full self-gravity. To maintain rea- 
sonable computational cost, numerical resolution in the rcj) 
plane is much reduced compared to razor-thin disc simula- 
tions (~ 4 — 6 cells per H compared to ~ 16 in 2D). Despite 
this, our plots in the rep plane closely resemble those ob- 
tained from 2D simulations. 

On the other hand, the low resolutions adopted here are 
unlikely to ca pture elliptic instabilities which may destroy 
vortices in 3D (|Lesur fe Papaloizoul 120091 . 1 20 id ). This might 
not be a serious issue when the condition for the vortex in- 
stability is maintained (by a planet in the present context). 
Also, the vortex grows on dynamical ti mescales whereas the 
ellipt ic instability takes much longer (|Lesur &; Papaloizoul 
120091 ). The initial development of vortices is not expected 
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to be suppressed by the elliptic instability (as found by 
iMeheut et aT1l2012bl ). 

The vortex instability produces s maller vortices with 
i ncrea sing self-gravity. According to iLesur &; Papaloizoul 
(2009), vortices with small aspect-ratios ( < 4) are strongly 
unstable in 3D, but note that their model is a local patch of 
a smooth disc without self-gravity, whereas we considered a 
gap edge in a global self- gravitating disc. 

Inclusion of self-gravity may change the stability prop- 
erties of vortices. In particular, we found tha t a vor- 
tex ca n flatten somewhat under its own weight. iLithwickl 
(2009) suggested vertical gravity helps to stabilise vor- 
tices in 3D. Vertical se/f-gravity can enhance this effect. 
ILithwickl found in local 3D simulations that 'tall' vor- 
tices are unstable whe reas 'short' vortices survive as in 2D 
(jGodon fe Livid Il999h . We may expect the more stratified 
vortices formed in self- gravitating disc s to be more stable 
than those in non- self- gravitating discs ([Barranco fc Marcus! 
120051). The elliptic instability is also weakened by stratifica- 
tion (|Lesur &; Papaloizoul [2009). so vertical self-gravity may 
also be stabilising in this respect. 

We have adopted the locally isothermal EOS for sim- 
plicity and direct comparison with previous 2D results. 
This EOS limits the 3D structure of densit y waves com- 
pared to thermally stratified discs (e 



iLubow fc Qgilvielll998l : lOdlvie &; Lubo 1 



.Lin et al.l Il990l : 
1999), which can 



cause refraction of waves out of the midplane. The lo- 
cally isoth ermal EOS represents the limit of efficient cooling 
(Boss |l998h . This might apply in optically thin regions of 
a dislQ but is violated if high densities d evelop, such as 
self- gravitating clumps ([Pickett et al.l l2000h . Our edge in- 
stabilities only reach moderate over-densities. Nevertheless, 
enhanced vertical stratification of the edge disturbances ob- 
served in our simulations are likely exaggerated by the EOS. 
Clearly, it is necessary to extend models of edge instabili- 
ties in self- gravitating discs to include an appropriate energy 
equation. 

Another important issue is vertical boundary condi- 
tions. We simulated a thin disc and imposed a reflecting up- 
per boundary to prevent mass loss from above. This setup 
may enhance the two-dimensionality of the problem. The 
vortex instability tends to involve the entire column of fluid, 
especially in the weak self- gravity regime. It is therefore a 
global instability in z. We suspect the spiral mode is less 
affected by vertical boundary conditions because the insta- 
bility tends to concentrate material at the midplane. Future 
simulations will consider varying the vertical domain size 
and upper disc boundary conditions. 

Our torque measurements indicate migration will be 
significantly affected by the vortex and spiral instabilities. 
Preliminary 3D simulations with a freely migrating planet 
have been performed. We recover the vortex-planet scat- 
terin g and the spiral - induced ou tward migration described 
bv lLin fe Papaloizoul <|2010L l2Q12h . The disc-planet torque is 
determined by the density field, and the above instabilities 
have density perturbations that either have weak vertical 
dependence or concentrated at the midplane. Thus, we be- 
lieve that at present, 2D simulations are more advantageous 



b For example, Cossi ns et al ] (120101 ) find optical depths r < 0.2 
beyond ~ 100AU in their models of protoplanetary discs. 



for studies focusing on migration, because high resolution is 
feasible and needed. 
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APPENDIX A: MODIFICATION TO VERTICAL 
STRUCTURE BY SELF-GRAVITY 

We describe a simple procedure to set up the vertical struc- 
ture of a locally isothermal, self- gravitating disc. We imagine 
setting up a non-self- gravitating disc, then slowly switch on 
the vertical force due to self-gravity. We expect the midplane 
density to increase at the expense of gas density higher in 
the atmosphere. It is assumed that the temperature profile 
remains unchanged. 

Vertical hydrostatic equilibrium between gas pressure, 
stellar gravity and self-gravity reads 



)\np 

dz 



dz dz 



(Al) 



Assuming a smooth radial density profile, we use the plane- 
parallel atmosphere approximation for the disc potential, 



d 2 $ 

dz 2 



4nGp. 



(A2) 



Next, we write the density field as 

p(R, z) = p N (R, z) x p(z] R) (A3) 

where pN is the density field corresponding to the non-self- 
gravitating disc: 
In p N _ 

~ ~~dz~' 



cL(R)- 



pN 



dz 
E(i2) 



(A4) 



■exp( 



-z 2 /2H 2 ) 



2ttH(R) 

= p N0 (R)e^p(-z 2 /2H 2 ), (A5) 

where pNo = Pn(R,z = 0) is the midplane density. The 
explicit expression for pN above assumes a thin disc. The 
function f3 describes the modification to the local density in 
order to be consistent with self- gravity. By construction, its 
governing equation is 

d 2 ln/3 



dz 2 



-4ttGpnP. 



Let 



X = m/3 - z 2 /2H 2 



1/2 



(A6) 

(A7) 
(A8) 



then the governing equation can be written in dimensionless 
form 

d 2 Y 

^| = -K-exp X , (A9) 
c 2 

K = ^2 . 

4tvGpnoH 2 

Note that K(R) is proportional to the local Keplerian 
Toomre parameter. Eq. IA9I can be further reduced to a 
first order differential equation, but this is unnecessary be- 
cause we pursue a numerical solution at the end. Appropri- 
ate boundary conditions are 

x (z = 0) = ln/5 0> (A10) 
Ox 



0. 



(All) 



flo(R) is the midplane density enhancement. To determine 
its value, we impose the surface density before and after 
modification by self- gravity to remain the same. Then we 
require 

F(M = \ — expx(6 A>)d6 - erf | 



= 0, (A12) 

where n is the number of scale-heights of the non-self- 
gravitating disc we originally considered. At a given cylin- 
drical radius R, we solve Eq. [Al2l using Newton-Raphson 
iteration. Each iteration involves integrating the governing 
ODE for x (Eq. IA9|) . At the end of the iteration, we have 
f3(z;R) and the midplane enhancement /3o(R). 

We comment that the procedure outlined above can be 
extended to polytropic discs. In this case, there is an addi- 
tional unknown — the new disc thickness after adjustment 
by self-gravity and an additional constraint — the density 
should vanish at the new disc surface. 



